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Abstract. We introduce a novel algorithm to perform event-driven simulations of hard rigid bodies of arbitrary shape, that 
relies on the evaluation of the geometric distance. In the case of a monodisperse system of uniaxial hard ellipsoids, we 
perform molecular dynamics simulations varying the aspect-ratio Xq and the packing fraction <j). We evaluate the translational 
Dtrans and the rotational D,-ot diffusion coefficient and the associated isodiffusivity lines in the ^ — Xq plane. We observe a 
decoupling of the translational and rotational dynamics which generates an almost perpendicular crossing of the Dtnms and 
Drot isodiffusivity lines. While the self intermediate scattering function exhibits stretched relaxation, i.e. glassy dynamics, 
only for large ^ and ~ 1, the second order orientational correlator C2{t) shows stretching only for large and small Xq 
values. We discuss these findings in the context of a possible pre-nematic order driven glass transition. 
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INTRODUCTION 

Particles interacting with only excluded volume interac- 
tion may exhibit a rich phase diagram, despite the ab- 
sence of any attraction. Simple non-spherical hard-core 
particles can form either crystalline or liquid crystalline 
ordered phases [? ], as first shown analytically by On- 
sager [? ] for rod-like particles. Although detailed phase 
diagrams of several hard-body (HB) shapes can be found 
in literature [????] less detailed information are avail- 
able about dynamics properties of hard-core bodies and 
their kinetically arrested states. 

The slowing down of the dynamics of the hard-sphere 
system on increasing the packing fraction ^ is well de- 
scribed by mode coupling theory (MCT) [? ], but on go- 
ing from spheres to non-spherical particles, non-trivial 
phenomena arise, due to the interplay between transla- 
tional and rotational degrees of freedom. The slowing 
down of the dynamics can indeed appear either in both 
translational and rotational properties or in just one of 
the two. Hard ellipsoids (HE) of revolution [? ? ] are one 
of the most prominent systems composed by hard body 
anisotropic particles. 

The equilibrium phase diagram, evaluated numerically 
two decades ago [? ], shows an isotropic fluid phase (I) 
and several ordered phases (plastic solid, solid, nematic 
N). The coexistence lines show a swallow-like depen- 
dence with a minimum at the spherical limit Xq = \ and 
a maximum at Xq k, 0.5 and Xqk,2 (cf. Figure [T]!. Ap- 
plication to HE [? ] of the molecular MCT (MMCT) [? ? 
] predicts also a swallow-like glass transition line. In ad- 
dition, the theory suggests that forZo ^ 0.5 and Xq ^ 2, 
the glass transition is driven by a precursor of nematic 



order, resulting in an orientational glass where the trans- 
lational density fluctuations are quasi-ergodic, except for 
very small wave vectors q. Within MCT, dynamic slow- 
ing down associated to a glass transition is driven by the 
amplitude of the static correlations. Since the approach 
of the nematic transition line is accompanied by an in- 
crease of the nematic order correlation function at ^ = 0, 
the non-linear feedback mechanism of MCT results in a 
glass transition, already before macroscopic nematic or- 
der occurs [? ]. In the arrested state, rotational motions 
become hindered. 

We perform event-driven (ED) molecular dynamics 
simulations, using a new algorithm [? ? ], that differently 
from previous ones [? ? ], relies on the evaluation of the 
distance between objects of arbitrary shape and that will 
be illustrated shortly. 

The outline of the manuscript is as follows. In the next 
section we illustrate shortly the new algorithm, that we 
proposed for simulating objects of arbitrary shape. Then 
in Sec. we illustrate the model used and we give all 
the details of the simulations we performed. In Sec. we 
show the results concerning the dynamics of the system 
investigated. The final section contains our conclusions. 

ALGORITHM DETAILS 

An event-driven algorithm for hard rigid 
bodies 

In an ED simulation the system is propagated until 
the next event occurs, where an event can be a collision 
between particles, a cell crossing (if linked lists are used). 



etc. All these events must be ordered with respect to time 
in a calendar and insertion, deletion and retrieving of 
events must be performed as efficiently as possible. 

One elegant approach has been introduced twenty 
years ago by Rapaport [? ], who proposed to arrange all 
events into an ordered binary tree (that is the calendar 
of events), so that insertion, deletion and retrieving of 
events can be done with an efficiency 0{logN), 0(1) and 
0{logN) respectively, where is the number of events in 
the calendar We adopt this solution to handle the events 
in our simulation; all the details of this method can be 
found in [? ]. Our ED algorithm can be schematized, as 
follows: 

1. Initialize the events calendar (predict collisions, 
cells crossings, etc.). 

2. Retrieve next event S" and set the simulation time to 
the time of this event. 

3. If the final time has been reached, terminate. 

4. If ^ is a collision between particles A and B then: 

(a) change angular and center-of-mass velocities 
of A and B (see [? ]). 

(b) remove from the calendar all the events (col- 
lisions, cell-crossings) in which A and B are 
involved. 

(c) predict and schedule all possible collisions for 
A and B. 

(d) predict and schedule the two cell crossings 
events for A and B. 

5. If is a cell crossing of a certain rigid body A: 

(a) update linked lists accordingly. 

(b) remove from calendar all events (collisions, 
cell-crossings) in which A is involved. 

(c) predict and schedule all possible collisions for 
A using the updated linked lists. 

6. go to step 2. 

All the details about linked lists can be found again in 
Ref. [? ], where an ED algorithm for hard spheres is illus- 
trated. We note also that in the case of HE, according to 
[? ], if the elongation of HE is big (i.e. one axes is much 
greater or much smaller than the others), the linked list 
method becomes inefficient at moderate and high densi- 
ties. In fact, given one ellipsoid A, using the linked lists, 
collision times of A with all ellipsoids in the same cell 
of A and with all the 26 adjacent cells have to be pre- 
dicted. In the case of rotationally symmetric ellipsoids, 
the number of time-of-collision predictions grows as Xq, 
if > 1 and as Xq, if < 1 [? ]. To overcome this 
problem, we developed also a new nearest neighbours list 
(NNL) method. At the begin of the simulation we build 
an oriented bounding parallelepiped (OBP) around each 
ellipsoid and we only predict collisions between ellip- 
soids having overlapping OBR In other words given an 



ellipsoid A, all the ellipsoids having overlapping OBP are 
the NNL of A. In addition the time-of-collisioni f, of each 
ellipsoid with its corresponding OBP is evaluated and the 
nnlQ of each HE is rebuilt at the time t = min,{f,}. The 
most time-consuming step in the case of hard rigid bod- 
ies is the prediction of collisions between HE, hence in 
the following a short description of the algorithm, that 
we developed, is given. 

Distance between two rigid bodies 

Our algorithm for predicting the collision time of two 
rigid bodies relies on the evaluation of the distance be- 
tween them. The surfaces of two rigid bodies, A and B, 
are implicitly defined by the following equations: 

/(x)=0 (la) 

^(x)=0 (lb) 

In passing we note that in the case of uniaxial hard ellip- 
soids, we have/(x,f) = '(x-ryi(f))XA(f)(x-ryi(f)) - 1 
and a similar expression holds for g(x), where r^i is the 
position of the center-of-mass of A and XA{t) is a matrix, 
that depends on the orientation of the A (see [? ]). It is 
assumed that, if a point x is inside the rigid body A (B), 
then /(x) < (g{x) < 0), while if it is outside /(x) > 
(g(x) > 0). Then the distance d between these two ob- 
jects can be defined as follows: 

d= min llx^ — xb|| (2) 

fM=o 

Equivalently, the distance between these two objects can 
be defined as the solution of the following set of equa- 
tions: 

fxA = -a^gxB (3a) 
/(xa)=0 (3b) 
^(xb)-O (3c) 

xa + PfxA = xb (3d) 

where x^ = (xa,3'a,Za), x^ = {xB,yB,ZB), Ua ^ and 

Sxb ~ '3^- '^l^- * l3b] i and ( l3c] i ensure that xa and x^ are 
points on A and B, eq.(l3al) provides that the normals to 
the surfaces are anti-parallel, and eq.(l3dl) ensures that the 
displacement of xa from xg is collinear to the normals 
of the surfaces. Equations ([3]) define extremal points of 
d; therefore, for two general HB these equations can 
have multiple solutions, while only the smallest one is 
the actual distance. To solve such equations iteratively is 



This is a particular case of the collision between rigid bodies. 



therefore necessary to start from a good initial guess of 
(x^,XB,a,/3) to avoid finding spurious solutions. 

In addition note that if two HB overlap slightly (i.e. the 
overlap volume is small) there is a solution with j3 < 
that is a measure of the inter-penetration of the two rigid 
bodies; we will refer to such a solution as the "negative 
distance" solution. 



2. Overestimating the rate of variation of the distance 
with respect to time (d{t)), refine the bracketing of 
the collision time obtained in 1 . 

3. Find the collision time to the best accuracy using a 
Newton-Raphson on a suitable set of equations for 
the contact point and the contact timeQ 



Newton-Raphson method for the distance 

The set of equations ^ can be conveniently solved 
by a Newton-Raphson (NR) method, as long as first and 
second derivatives of /(x^) and ^(xb) are well defined. 
This method, provided that a good initial guess has been 
supplied, very quickly reaches the solution because of its 
quadratic convergence [? ]. If we define: 



F(y) = 



Eqs. (O become: 



F(y) = 



(4) 



(5) 



wherey= (xA,XB,a,j3). 

Given an initial point yo we build a sequence of points 
converging to the solution as follows: 



y,-+i=y, + r'F(yO 



(6) 



where J is the Jacobian of F,i.e. J 

The matrix inversion, required to evaluate J^', can be 
done making use of a standard LU decomposition [? ]. 
This decomposition is of order /3, where is the 
number of equations (8 in the present case). Finally we 
note that the set of 8 equations ^ can be also reduced to 
5 equations, ehminating x^ or x^, using Eq. (l3d] i . 



Prediction of the time-of-coUision 

The collision (or contact) time of two rigid bodies A 
and B is, the smallest time tc, such that d{tc) = 0, where 
d{t) is the distance as a function of time between A and 
B. For finding tc we perform the following steps: 

1 . Bracket the contact time using centroids (this tech- 
nique will not be discussed here, see [? ? ] for the 
detailsfl 



Bracketing of the contact time overestimating d{t). 

It can be proved that an overestimate of the rate of 
variation of the distance is the following: 



d{t) < ||va - vbII + ||wa||La + IIwbIILb 



(7) 



where the dot indicates the derivation with respect to 
time, Ta and are the centers of mass of the two rigid 
bodies, va and \b are the velocities of the centers of 
mass, and the lengths L^,, Lb are such that 



La > max {||r'-rA||} 

/(r')<0 

Lb > max {||r'-rB||} 



<o 



(8a) 



(8b) 



Using this overestimate of d{t), that will be called 
dmax from now on, an efficient strategy, to refine the 
bracketing of the contact time, is the following: 

1 . If f 1 and f2 bracket the solution, set t —ti. 

2. Evaluate the distance d{t) at time t. 

3. Choose a time increment Af as follows: 



At 



d(t) 



if d{t) > Eci; 
otherwise. 



(9) 



where <C min{LA,LB} 
Evaluate the distance at time t + At. 
If d{t + At) < and d{t) > 0, then fi = t and t2 = 
t + At, find the collision time/point via NR (see Sec. 
) and terminate (collision will occur after fi). 
if both < d{t + At) < £,/ and < d{t) < £,/, there 
could be a "grazing" collision between t and t + At 
[? ] (distance is first decreasing and then increas- 
ing). To look for a possible collision, evaluate the 
distance d{t + At /2) and perform a quadratic inter- 
polation of the 3 points ( {t,d{t)), {t +At/2,d{t + 
At/2), {t +At,d{t +At) ). If the resulting parabola 
has zeros, set t\ to the smallest zero (first collision 
will occur near the smallest zero), and find the col- 
lision time/point via NR and terminate (see Sec. ). 



- A centroid of a given ellipsoid A with center is the smallest sphere 
centered at r^^ that that encloses A 



' Alternatively you can use any one-dimensional root-finder for the 
equation d{t) = 0. 



7. Increment time by f ^ f + Af. 

8. if f > f2 terminate (no collision has been found). 

9. Go to step 2. 

Finally we note that if the quadratic interpolation fails, 
the collision will be missed, anyway if e^/ is enough small 
all "grazing" collisions will be properly handled, i.e. all 
collisions will be correctly predicted. 



Set of equations to find the contact time and the contact 
point 

The contact time, to the best possible accuracy, can be 
found solving the following equations: 



/x(x,f) = -a2gx(x,f) 



/(x,f) 








(10a) 

(10b) 
(10c) 



where now / and g depend also on time because the 
two objects move, and the independent variables are the 
contact time and the contact point. Again, a good way 
to solve such a system is using the NR method, very 
similarly to what we did for evaluating the geometric 
distance. NR for Eqs. (fTol l is again very unstable unless a 
very good initial guess is provided, but the bracketing 
evaluated using dmax is sufficiently accurate, provided 
that is enough small (typically ^ l/s — 4 is a good 
choice to give a good initial guess for this NR and to 
avoid "grazing" collisions). 

Finally we want to stress that this algorithm will be 
exploited fully, when hard bodies of arbitrary shape will 
be simulated, and that with minor changes it can also 
be used to simulate hard bodies decorated with attractive 
spots, as it has been done in the past for the specific case 
of hard-spheres [? ? ]. Work along these directions is on 
the way, and in particular a system of super-ellipsoids 
has been already successfully simulated using the present 
algorithm. 



METHODS 

We perform an extended study of the dynamics of 
monodisperse HE in a wide window of and Xq val- 
ues, extending the range of Xq previously studied [? ]. 
We specifically focus on establishing the trends leading 
to dynamic slowing down in both translations and rota- 
tions, by evaluating the loci of constant translational and 
rotational diffusion. These lines, in the limit of vanishing 
diffusivities, approach the glass-transition lines. We also 
study translational and rotational correlation functions, 
to search for the onset of slowing down and stretching 



in the decay of the correlation. We simulate a system of 
= 512 ellipsoids at various volumes V = in a cu- 
bic box of edge L with periodic boundary conditions. We 
chose the geometric mean of the axis / — ifal? as unit 
of distance, the mass m of the particle as unit of mass 
{m — 1) and hgT — 1 (where kg is the Boltzmann con- 
stant and T is the temperature) and hence the correspond- 
ing unit of time is \fmP^Jk^ . The inertia tensor is cho- 
sen as Ix — ly = 2mr^/5, where r — min{a,b}. The value 
of the L component is irrelevant [? ], since the angular 
velocity along the symmetry (z-) axis of the HE is con- 
served. We simulate a grid of more than 500 state points 
at different and ^.To create the starting configuration 
at a desired (j), we generate a random distribution of ellip- 
soids at very low (p and then we progressively decrease 
L up to the desired (j). We then equilibrate the configu- 
ration by propagating the trajectory for times such that 
both angular and translational correlation functions have 
decayed to zero. Finally, we perform a production run at 
least 30 times longer than the time needed to equilibrate. 
For the points close to the I-N transition we check the 
nematic order by evaluating the largest eigenvalue S of 
the order tensor Q [? ], whose components are: 



3 1 1 



(11) 



where aj3 € {xjjjz}, and the unit vector {ui{t))a is the 
component a of the orientation (i.e. the symmetry axis) 
of ellipsoid / at time t. The largest eigenvalue S is non- 
zero if the system is nematic and if it is isotropic. In 
the following, we choose the value S = 0.3 as criteria to 
separate isotropic from nematic states. 



RESULTS AND DISCUSSION 
Isodiffusivity lines 

From the grid of simulated state points we build a cor- 
responding grid of translational (Aran*) and diffusional 
(Drut) coefficients, defined as: 



A™„,= iimlEMO^i« 

t^+^N^ 6t 



Dr. 



lim — 

r^+~ A^ ' 



lA*,: 



4f 



(12) 



(13) 



where AO,- = Jq (Oidt, x, is the position of the center 
of mass and O; is the angular velocity of ellipsoid /. 
By proper interpolation, we evaluate the isodiffusivity 
lines, shown in Fig. \T\ Results show a striking decou- 
pling of the translational and rotational dynamics. While 



the translational isodiffusivity lines mimic the swallow- 
like shape of the coexistence between the isotropic liq- 
uid and the crystalline phases (as well as the MMCT 
prediction for the glass transition [? ]), rotational isod- 
iffusivity lines reproduce qualitatively the shape of the 
I-N coexistence. As a consequence of the the swallow- 
like shape, at large fixed 0, Dirans increases by increas- 
ing the particle's anisotropy, reaching its maximum at 
X() « 0.5 and Xq « 2. Further increase of the anisotropy 
results in a decrease of Dtrans- For all Xq, an increase 
of ^ at constant Xq leads to a significant suppression of 
Dtrans, demonstrating that Dtrans is controlled by pack- 
ing. The iso-rotational lines are instead mostly controlled 
by Xq, showing a progressive slowing down of the ro- 
tational dynamics independently from the translational 
behavior This suggests that on moving along a path of 
constant Dtrans, it is possible to progressively decrease 
the rotational dynamics, up to the point where rotational 
diffusion arrest and all rotational motions become hin- 
dered. Unfortunately, in the case of monodisperse HE, 
a nematic transition intervenes well before this point is 
reached. It is thus stimulating to think about the possibil- 
ity of designing a system of hard particles in which the 
nematic transition is inhibited by a proper choice of the 
disorder in the particle's shape/elongations. We note that 
the slowing down of the rotational dynamics is consis- 
tent with MMCT predictions of a nematic glass for large 
Xq he [? ], in which orientational degrees of freedom 
start to freeze approaching the isotropic-nematic transi- 
tion line, while translational degrees of freedom mostly 
remain ergodic. 



Orientational correlation function 

To support the possibility that the slowing down of 
the dynamics on approaching the nematic phase orig- 
inates from a close-by glass transition, we evaluate 
the self part of the intermediate scattering function 
Fseif{q,t) = l(£^-e'i(''j(')-''j(0))) and the second order 
orientational correlation function €2(1) defined as [? ] 
Ciit) = (P2(cos0(f))), where P2{x) = {3x^-l)/2 and 
0(f) is the angle between the symmetry axis at time t 
and at time 0. The €2(1) rotational isochrones are found 
to be very similar to rotational isodiffusivity lines. 

These two correlation functions never show a clear 
two-step relaxation decay in the entire studied region, 
even where the isotropic phase is metastable, since the 
system can not be significantly over-compressed. As for 
the well known hard-sphere case, the amount of over- 
compressing achievable in a monodisperse system is 
rather limited. This notwithstanding, a comparison of the 
rotational and translational correlation functions reveals 
that the onset of dynamic slowing down and glassy dy- 




FIGURE 1. [Color online] Isodiffusivity lines. Solid lines 
are isodiffusivity lines from translational diffusion coefficients 
Dtrans and dashed lines are isodiffusivities lines from rotational 
diffusion coefficients Dror- Arrows indicate decreasing diffu- 
sivities. Left and right arrows refer to rotational diffusion co- 
efficients. Diffusivities along left arrow are: 1.5, 0.75, 0.45, 
0.3, 0.15. Diffusivities along right arrow are: 1.5, 0.75, 0.45, 
0.3, 0.15, 0.075, 0.045. Central arrow refers to translational 
diffusion coefficients, whose values are: 0.5, 0.3, 0.2, 0.1, 0.04, 
0.02. Thick Long-dashed curves are coexistence curves of all 
first order phase transitions in the phase diagram of HE eval- 
uated by Frenkel and Mulder (FM) [? ] Solid lines are coex- 
istence curves for the I-N transition of oblate and prolate el- 
lipsoids, obtained analytically by Tijpto-Margo and Evans [? ] 
(TME). 



namics can be detected by the appearance of stretching. 

We note that F^eif shows an exponential behaviour 
close to the I-N transition (Xq = 3.2, 0.3448) on the 
prolate and oblate side, in agreement with the fact that 
translational isodiffusivities lines do not exhibit any pe- 
culiar behaviour close to the I-N line [? ]. Only when 

« 1, Fggif develops a small stretching, consistent with 
the minimum of the swallow-like curve observed in the 
fluid-crystal line [? ? ], in the jamming locus as well as 
in the predicted behavior of the glass line for HE [? ] 
and for small elongation dumbbells [? ? ]. Opposite be- 
havior is seen for the case of the orientational correla- 
tors. C2 shows stretching at large anisotropy, i.e. at small 
and large Xq values, but decays within the microscopic 
time for almost spherical particles. In this quasi-spherical 
limit, the decay is well represented by the decay of a free 
rotator [? ? ]. Previous studies of the rotational dynamics 
of HE [? ] did not report stretching in C2, probably due 
to the smaller values of Xq previously investigated and to 
the present increased statistic which allows us to follow 
the full decay of the correlation functions. 

In summary C2 becomes stretched approaching the I-N 
transition while F^^i f remains exponential on approach- 
ing the transition. To quantify the amount of stretching in 
C2, we fit it to the function A exp [-(?/%,) '^''2] (stretched 
exponential) for several state points and we show in Fig. 
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FIGURE 2. [Color online] jSc, and Tc, are obtained from fits 
of C2 to a stretched exponential for <j) = 0.40,0.45 and 0.50. 
Top: Tc, as a function of Xq. Bottom: /3c, as a function of Xq. 
The time window used for the fits is chosen in such a way 
to exclude the microscopic short times ballistic relaxation. For 
0.588 <Xq < 1.7 the orientational relaxation is exponential. 



12] the Xq dependence of and for three different 
values of 0. In all cases, slowing down of the character- 
istic time and stretching increases progressively on ap- 
proaching the I-N transition. 



CONCLUSIONS 

In summary we have applied a novel algorithm for simu- 
lating hard bodies to investigate the dynamics properties 
of a system of monodisperse HE and we have shown that 
clear precursors of dynamic slowing down and stretch- 
ing can be observed in the region of the phase diagram 
where a (meta)stable isotropic phase can be studied. De- 
spite the monodisperse character of the present system 
prevents the possibility of observing a clear glassy dy- 
namics, our data suggest that a slowing down in the ori- 
entational degrees of freedom — driven by the elonga- 
tion of the particles — is in action. The main effect of 
this shape-dependent slowing down is a decoupling of 
the translational and rotational dynamics which gener- 
ates an almost perpendicular crossing of the Dtrans and 
Diot isodiffusivity lines. This behavior is in accordance 
with MMCT predictions, suggesting two glass transition 
mechanisms, related respectively to cage effect (active 
for 0.5 ^Xq ^2) and to pre-nematic order {Xq ^ 0.5, 
^o^2)T?]. ^ 
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